##Replication code for "disentangling the combined effects of PEC and legislative polarisation"
library(ggplot2)
library(ggtext)
library(gridExtra)
library(haven)
library(marginaleffects)
library(mclogit)
library(patchwork)
library(ragg)
library(readxl)
library(rmapshaper)
library(rnaturalearth)
library(rnaturalearthdata)
library(showtext)
library(sf)
library(stargazer)
library(survival)
library(tidyverse)

masters <- read.csv("govfor_rr.csv", sep = ",") %>%
           filter(presidential_election == 1 | parliamentary_election == 1) %>%
           mutate(pec = as.factor(pec)) 


##Dataviz Fig1----
#Data-Wrangling
masters_last_election <- masters %>%
                         group_by(country_name_long, start_date, election_id, polariz_st) %>%
                         summarise(n = n()) %>%
                         ungroup() %>%
                         filter(election_id %in% c("101019", "102009", "103009", 
                                                   "104008", "105012", "106017",
                                                   "109010", "110008", "111015",
                                                   "112008", "114005", "107011")) %>%
                         filter(country_name_long != "Bolivia") %>%
                         group_by(country_name_long, election_id, polariz_st) %>%
                         summarise(n = n()) %>%
                         ungroup() %>%
                         select(-n) %>%
                         rename(name_long = country_name_long)


masters_average <- masters %>%
                    group_by(country_name_long, election_id, polariz_st) %>%
                    summarise(n = n()) %>%
                    ungroup() %>%
                    select(-n) %>%
                    group_by(country_name_long) %>%
                    summarise(mean_polariz = mean(polariz_st, na.rm = T)) %>%
                    rename(name_long = country_name_long)

world <- ne_countries(scale = "medium", 
                      returnclass = "sf")

nowayoutofhere <- world %>%
                  filter(region_wb == "Latin America & Caribbean") 


nowayoutofhere_points <- st_centroid(nowayoutofhere) %>%
                         cbind(st_coordinates(st_centroid(nowayoutofhere$geometry)))

pneuma <- nowayoutofhere_points %>%
          st_drop_geometry() %>%
          select(name_long, X, Y) %>%
          left_join(masters_average) %>%
          left_join(masters_last_election) %>%
          left_join(nowayoutofhere) 

#Fonts
#You should have the fonts below to reproduce all visualisations of this script.
#They are available here: https://fonts.google.com/specimen/Barlow+Semi+Condensed
font_add(family = "regular", "Barlow Semi Condensed-Regular.ttf")
font_add(family= "bold", "BarlowSemiCondensed-SemiBold.ttf")
showtext_auto()   

#Plots
average <- ggplot() +
  geom_sf(data = pneuma, aes(geometry = geometry), inherit.aes = FALSE) +
  geom_sf(data = pneuma %>% filter(!is.na(mean_polariz)), 
          aes(fill = mean_polariz, geometry = geometry), 
          size = .15) +
  geom_sf(data =  pneuma %>% filter(is.na(mean_polariz)) %>% 
            mutate(mean_polariz = ifelse(is.na(mean_polariz), 0, mean_polariz)),
          aes(fill = mean_polariz, geometry = geometry),
          size = .15, 
          linetype = 2) +
  coord_sf(expand = FALSE) +
  scale_fill_gradientn(breaks = c(0, 1.5, 2.5, 3.5), 
                       colours = c("#ffffff", "#dbdbdd", "#c4c6c8", "#a7a9ab", "#606163", "#333333"), 
                       labels = c("NA", 1.5, 2.5, 3.5)
  ) + 
  xlab("") + ylab("") +
  labs(title = "Average Legislative Polarisation in Latin America") +
  theme(##Background
    panel.border=element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "#ffffff", linetype = 'blank'),
    panel.background = element_rect(fill = "#ffffff"),
    ##Axes
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    ##Legend
    legend.text = element_text(size=15, family="regular"),
    legend.position = c(0.1, 0.2),
    legend.key.size = unit(1.0, "cm"),
    #legend.key = element_rect(fill = "white", colour = "#323233", linewidth = 0.5),
    legend.spacing.y = unit(0, "cm"),
    ##Titles
    plot.title=element_markdown(family="bold", hjust=0.2, vjust = 0.0, size=18, color="black"),
    plot.title.position = "plot",
  ) +
  guides(fill = guide_legend(title = "", 
                             reverse = T,
                             byrow = TRUE,
                             override.aes = list(linetype = c("solid",
                                                              "solid",
                                                              "solid",
                                                              "dashed"))),
         list(legend.key = element_rect(colour = c("black",
                                                   "black",
                                                   "black",
                                                   "white")))) 

last <- ggplot() +
  geom_sf(data = pneuma, aes(geometry = geometry), inherit.aes = FALSE) +
  geom_sf(data = pneuma %>% filter(!is.na(polariz_st)), 
          aes(fill = polariz_st, geometry = geometry), 
          size = .15) +
  geom_sf(data =  pneuma %>% filter(is.na(polariz_st)) %>% 
            mutate(polariz_st = ifelse(is.na(polariz_st), 0, polariz_st)),
          aes(fill = polariz_st, geometry = geometry),
          size = .15, 
          linetype = 2) +
  coord_sf(expand = FALSE) +
  scale_fill_gradientn(breaks = c(0, 1, 2, 3, 4), 
                       colours = c("#ffffff", "#dbdbdd", "#c4c6c8", "#a7a9ab", "#606163", "black"), 
                       labels = c("NA", 1, 2, 3, 4)
  ) + 
  xlab("") + ylab("") +
  labs(title = "Legislative Polarisation in Latin America from 2016 to 2019") +
  theme(##Background
    panel.border=element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "#ffffff", linetype = 'blank'),
    panel.background = element_rect(fill = "#ffffff"),
    ##Axes
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    ##Legend
    legend.text = element_text(size=15, family="regular"),
    legend.position = c(0.1, 0.2),
    legend.key.size = unit(1.0, "cm"),
    #legend.key = element_rect(fill = "white", colour = "black", linewidth = 0.2),
    legend.spacing.y = unit(0, "cm"),
    ##Titles
    plot.title=element_markdown(family="bold", hjust=0.2, vjust = 0.0, size=18, color="black"),
    plot.title.position = "plot",
  ) +
  guides(fill = guide_legend(title = "", 
                             reverse = T,
                             byrow = TRUE,
                             override.aes = list(linetype = c("solid",
                                                              "solid",
                                                              "solid",
                                                              "solid",
                                                              "dashed"))))

grid.arrange(average, last, ncol=2)
##Dataviz Fig2----
#Data-Wrangling
#Before making this plot, you should download the V-Party data from the v-dem repository: https://www.v-dem.net/data/v-party-dataset/
vdem <- read.csv("V-Dem-CPD-Party-V2.csv", sep=",") %>%
        filter(country_name == "Chile") %>%
        filter(year == 2009) %>%
        select(v2paenname, v2pashname, v2paseatshare, v2pariglef_osp) %>%
        mutate(ideology = 1+(19/6)*v2pariglef_osp)

pinera <- vdem %>%
           drop_na() %>%
           mutate(ideology_vd2 = ideology/2) %>%
           filter(v2pashname == "RN" | v2pashname == "UDI") %>%
           mutate(v2paseatshare = v2paseatshare/100)

pinera_all  <- vdem %>%
               drop_na() %>%
               mutate(ideology_vd2 = ideology/2) %>%
               mutate(v2paseatshare = v2paseatshare/100)

#Theme
mytheme2 <- theme(
            ##Titles
            plot.title=element_text(family="bold", hjust=0.6, vjust = 0.0, size=17, color="black"),
            plot.title.position = "plot",
            plot.subtitle = element_markdown(family="bold", size=20, hjust=0.5, color="black"),
            plot.caption= element_text(family="bold", size=15, color="black", hjust=0.5),
            plot.caption.position = "plot",
            ##Background
            panel.border=element_blank(),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            plot.background = element_rect(fill = "transparent", linetype = 'blank'),
            panel.background = element_rect(fill = "transparent"),
            #plot.margin=ggplot2::margin(0.5, 0.5, 0.5, 0.5, "in"),
            #Axes
            axis.ticks.length=unit(0.15, "cm"),
            axis.ticks = element_blank(),
            axis.line = element_blank(),
            axis.title = element_text(size=15, family="regular", color="black"),
            axis.text.y = element_text(size=15, family="regular", color="black"),
            axis.text.x = element_text(size=15, family="regular", color="black"),
            #no legend
            legend.position = "none")

a <- pinera  %>% ggplot(aes(ideology_vd2, v2paseatshare, label = v2pashname)) +
                 geom_point(aes(colour = "black", fill = "black"), size = 6, alpha = 0.85) +
                 scale_color_manual(values = c("black", "black")) +
                 scale_fill_manual(values = c("black", "black")) +
                 geom_text(aes(x=ideology_vd2, y=v2paseatshare+0.03, label = v2pashname), size = 4) +
                 coord_cartesian(expand=T) +
                 scale_x_continuous(limits = c(0, 11),
                                    breaks = seq(1, 10, by = 1)) +
                 scale_y_continuous(breaks = seq(00.10, 00.50, by = 00.10),
                                    limits = c(0.00, 00.50),
                                    labels = scales::percent_format(accuracy = 1)) +
                 ylab("% of seats in the Lower Chamber") + xlab("Left-Right Dimension") +
                 labs(title = "Piñera's PEC and First Cabinet") +
                 mytheme2

b <- pinera_all %>% ggplot(aes(ideology_vd2, v2paseatshare, label = v2pashname)) +
                     geom_point(aes(colour = "black", fill = "black"), size = 6, alpha = 0.85) +
                     scale_color_manual(values = c("black", "black")) +
                     scale_fill_manual(values = c("black", "black")) +
                     geom_text(aes(x=ideology_vd2, y=v2paseatshare+0.03, label = ifelse(v2pashname %in% c("PS", "PDC", "RN", "UDI"), v2pashname, "")), size = 4) +
                     geom_text(aes(x=ideology_vd2-0.8, y=v2paseatshare+0.03, label = ifelse(v2pashname=="PPD", v2pashname, ""), size = 4)) +
                     geom_text(aes(x=ideology_vd2+0.5, y=v2paseatshare+0.03, label = ifelse(v2pashname=="PRSD", v2pashname, ""), size = 4)) +
                     coord_cartesian(expand=T) +
                     scale_x_continuous(limits = c(0, 11),
                                        breaks = seq(1, 10, by = 1)) +
                     scale_y_continuous(breaks = seq(00.10, 00.50, by = 00.10),
                                        limits = c(0.00, 00.50),
                                        labels = scales::percent_format(accuracy = 1)) +
                     ylab("% of seats in the Lower Chamber") + xlab("Left-Right Dimension") +
                     labs(title = "Chilean Party System as of 2010") + 
                     mytheme2

grid.arrange(a, b, ncol=2)


##Dataviz Fig3----
#Before making this plot, you should download the dataset used for Borges et al. (forthcoming) "The Recasting of the Latin American 
#Right: Polarization and Conservative Reactions". At the time of writing, the book and the dataset are not publicly available. However, they should be available as soon as they
#are officialy released. Have a check at https://andreborges.org/datasets/.

`%notin%` <- Negate(`%in%`)

nr_s <- read.csv("NewRight_partylevel_dataset_v1.2.csv", sep = ";")
obs <- nr_s %>%
  filter(countrycode == "Brazil2010") %>%
  select(party, legseat, ideo_BG) %>%
  drop_na() %>%
  mutate(ideology = ideo_BG/2)

polas <- obs %>%
  mutate(all = sum(legseat)) %>%
  mutate(pct = legseat/all) %>%
  mutate(mlr = pct*ideology) %>%
  mutate(general = sum(mlr)) %>%
  mutate(polar = (ideology-general)/9.5 * 10) %>%
  mutate(polar = sqrt((polar*pct)^2)) %>%
  mutate(polar = sum(polar)) %>%
  select(polar) %>% unique() %>% pull()

rousseff <- obs %>%
  filter(party %in% c("PCdoB", "PT", "PDT", "PMDB", "PR", "PRB",
                      "PSB", "PSC")) %>%
  mutate(legseat = legseat/100)

rousseff_all  <- obs %>%
  mutate(legseat = legseat/100)

#Theme
mytheme2 <- theme(
  ##Titles
  plot.title=element_text(family="bold", hjust=0.6, vjust = 0.0, size=17, color="black"),
  plot.title.position = "plot",
  plot.subtitle = element_markdown(family="bold", size=20, hjust=0.5, color="black"),
  plot.caption= element_text(family="bold", size=15, color="black", hjust=0.5),
  plot.caption.position = "plot",
  ##Background
  panel.border=element_blank(),
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  plot.background = element_rect(fill = "transparent", linetype = 'blank'),
  panel.background = element_rect(fill = "transparent"),
  #plot.margin=ggplot2::margin(0.5, 0.5, 0.5, 0.5, "in"),
  #Axes
  axis.ticks.length=unit(0.15, "cm"),
  axis.ticks = element_blank(),
  axis.line = element_blank(),
  axis.title = element_text(size=15, family="regular", color="black"),
  axis.text.y = element_text(size=15, family="regular", color="black"),
  axis.text.x = element_text(size=15, family="regular", color="black"),
  #no legend
  legend.position = "none")

a <- rousseff  %>% mutate(ministerial_post = ifelse(party %notin% c("PRB", "PSC"), 1, 0)) %>%
                   ggplot(aes(ideology, legseat, label = party, colour = as.factor(ministerial_post), fill = as.factor(ministerial_post))) +
                   geom_point(size = 6, alpha = 0.85) +
                   scale_color_manual(values = c("#808080", "black")) +
                   scale_fill_manual(values = c("#808080", "black")) +
                   geom_text(aes(x=ideology, y=legseat+0.03, label = ifelse(party %in% c("PCdoB", "PMDB", "PR", "PRB", "PSB",
                                                                                         "PSC", "PT"), party, "")), size = 4, family = "regular") +
                   geom_text(aes(x=ideology+0.55, y=legseat+0.03, label = ifelse(party %in% c("PDT"), party, "")), size = 4, family = "regular") +
                   coord_cartesian(expand=T) +
                   scale_x_continuous(limits = c(0, 11),
                                      breaks = seq(1, 10, by = 1)) +
                   scale_y_continuous(breaks = seq(00.10, 00.50, by = 00.10),
                                      limits = c(-0.02, 00.50),
                                      labels = scales::percent_format(accuracy = 1)) +
                   ylab("% of seats in the Lower Chamber") + xlab("Left-Right Dimension") +
                   labs(title = "Rousseff's Pre-Electoral Coalition") +
                   mytheme2

b <- rousseff_all %>% mutate(ministerial_post = ifelse(party %in% c("PP"), 0, 1)) %>%
                      ggplot(aes(ideology, legseat, label = party, colour = as.factor(ministerial_post), fill = as.factor(ministerial_post))) +
                      geom_point(size = 6, alpha = 0.85) +
                      scale_color_manual(values = c("#a9a9a9", "black")) +
                      scale_fill_manual(values = c("#a9a9a9", "black")) +
                      geom_text(aes(x=ideology, y=legseat+0.03, label = ifelse(party %in% c("PT", "PMDB", "PSDB",
                                                                                            "PCdoB", "PSB", "PSC", "PR"), party, "")), size = 4,
                                family = "regular") +
                      geom_text(aes(x=ideology, y=legseat-0.03, label = ifelse(party == "PRB", party, ""), size = 4, family = "regular")) +
                      geom_text(aes(x=ideology-0.05, y=legseat-0.03, label = ifelse(party == "PV", party, ""), size = 4, family = "regular")) +
                      geom_text(aes(x=ideology+0.55, y=legseat+0.03, label = ifelse(party %in% c("PDT"), party, "")), size = 4, family = "regular") +
                      geom_text(aes(x=ideology+0.55, y=legseat+0.03, label = ifelse(party %in% c("PFL"), party, "")), size = 4, family = "regular") +
                      geom_text(aes(x=ideology+0.3, y=legseat+0.03, label = ifelse(party %in% c("PTB"), party, "")), size = 4, family = "regular") +
                      geom_text(aes(x=ideology, y=legseat-0.03, label = ifelse(party %in% c("PP"), party, "")), size = 4, family = "regular") +
                      geom_text(aes(x=ideology, y=legseat-0.03, label = ifelse(party %in% c("PPS"), party, "")), size = 4, family = "regular") +
                      coord_cartesian(expand=T) +
                      scale_x_continuous(limits = c(0, 11),
                                         breaks = seq(1, 10, by = 1)) +
                      scale_y_continuous(breaks = seq(00.10, 00.50, by = 00.10),
                                         limits = c(-0.02, 00.50),
                                         labels = scales::percent_format(accuracy = 1)) +
                      ylab("% of seats in the Lower Chamber") + xlab("Left-Right Dimension") +
                      labs(title = "Brazilian Party System as of 2010") + 
                      mytheme2
  
  
  
  

grid.arrange(a, b, ncol=2)
##Dataset Summary----
des_stats <- masters %>%
             filter(empgov == 1) %>%
             drop_na(polariz_st) %>%
             select(country_name_long, cab_id, year) %>% 
             unique() %>%
             group_by(country_name_long) %>%
             summarise(n = n()) %>%
             ungroup() %>%
             mutate(total = sum(n)) %>%
             mutate(percentage = (n/total)*100)

##Base Models----
freude <- masters %>%
          filter(majsit == 0) 

#interaction
model_interact <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=freude
)

#all_presidents
model_interact_majority <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters
)

stargazer(model_interact, model_interact_majority, 
          type = "latex",
          covariate.labels = c("Minority", "Number of Parties", 
                               "Ideological Division", "Median Party", 
                               "Extreme Parties", "Runner-up Party",
                               "Pre-Electoral Coalition (PEC)", " PEC * Legislative Polarisation"),
          column.labels   = c("Minority Presidents", "All Presidents"),
          column.separate = c(1, 1),
          dep.var.labels.include = F,
          #These lines come from re-running the models with clogit from the package Survival. I include them here to make them comparable to the results
          #reported by Indridason (2011) and Freudenreich (2016). All the following regression outputs include similar lines. 
          add.lines=list(c("Number of Alternative Cabinets","147,736", "149,452"),
                         c("Log Likelihood", "-251.904", "-302.924")))

##Marginal Effects----
dt <- plot_cme(model_interact, variables = "pec", condition = "polariz_st", type = "link", draw = F)

plot_cme(model_interact, variables = "pec", condition = "polariz_st", type = "link") +
  ggplot2::geom_line(data = dt, ggplot2::aes(x = polariz_st, 
                                             y = estimate), size = 0.8) +
  geom_rug(aes(x = polariz_st), data = freude) +
  ggplot2::geom_ribbon(data = dt, alpha = 0.3, ggplot2::aes(x = polariz_st, 
                                                            y = estimate, 
                                                            ymin = conf.low, 
                                                            ymax = conf.high)) +
  
  geom_hline(yintercept = 0) +
  scale_y_continuous(breaks = seq(-5, 15, by =5))+
  xlab("Legislative Polarisation") + ylab("Marginal Effect of Pre-Electoral Coalitions on Cabinet Formation") +
  theme(
    ##Background
    panel.border = element_rect(color = "black", fill = NA),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "#ffffff", linetype = 'blank'),
    panel.background = element_rect(fill = "#ffffff"),
    ##Axes
    axis.title = element_text(size=15, family="regular", color="black"),
    axis.text.y = element_text(size=10, family="regular", color="black"),
    axis.text.x = element_text(size=10, family="regular", color="black")) -> du


et <- plot_cme(model_interact, variables = "polariz_st", condition = "pec", type = "link", draw = F)
plot_cme(model_interact, variables = "polariz_st", condition = "pec", type = "link") +
  xlab("Pre-Electoral Coalition") + ylab("Marginal Effect of Legislative Polarisation on Cabinet Formation") +
  theme(
    ##Background
    panel.border = element_rect(color = "black", fill = NA),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "#ffffff", linetype = 'blank'),
    panel.background = element_rect(fill = "#ffffff"),
    ##Axes
    axis.title = element_text(size=15, family="regular", color="black"),
    axis.text.y = element_text(size=10, family="regular", color="black"),
    axis.text.x = element_text(size=10, family="regular", color="black")) -> dg

dg + du

##Supplementary Material----
#Herrera's example
#Once again, to reproduce the visualisation below, you must have access to the dataset used for Borges at al. (forthcoming)
nr_venezuela <- nr_s %>%
                filter(countrycode == "Venezuela1978") %>%
                select(party, legseat, ideo_BG) %>%
                mutate(ideo_BG = ifelse(party == "URD", 14.80, ideo_BG)) %>%
                drop_na() %>%
                filter(legseat >= 1.0) %>%
                mutate(ideology = ideo_BG/2)

polas <- nr_venezuela %>%
         mutate(all = sum(legseat)) %>%
         mutate(pct = legseat/all) %>%
         mutate(mlr = pct*ideology) %>%
         mutate(general = sum(mlr)) %>%
         mutate(polar = (ideology-general)/9.5 * 10) %>%
         mutate(polar = sqrt((polar*pct)^2)) %>%
         mutate(polar = sum(polar)) %>%
         select(polar) %>% unique() %>% pull()

herrera <- nr_venezuela %>%
           filter(party %in% c("COPEI", "URD")) %>%
           mutate(legseat = legseat/100)

herrera_all  <- nr_venezuela %>%
  mutate(legseat = legseat/100)


a <- herrera  %>% mutate(ministerial_post = ifelse(party %notin% c("URD"), 1, 0)) %>%
                  ggplot(aes(ideology, legseat, label = party, colour = as.factor(ministerial_post), fill = as.factor(ministerial_post))) +
                  geom_point(size = 6, alpha = 0.85) +
                  scale_color_manual(values = c("#808080", "black")) +
                  scale_fill_manual(values = c("#808080", "black")) +
                  geom_text(aes(x=ideology, y=legseat+0.03, label = party, size = 4, family = "regular"))+
                  coord_cartesian(expand=T) +
                  scale_x_continuous(limits = c(0, 11),
                                     breaks = seq(1, 10, by = 1)) +
                  scale_y_continuous(breaks = seq(00.10, 00.50, by = 00.10),
                                     limits = c(-0.02, 00.50),
                                     labels = scales::percent_format(accuracy = 1)) +
                  ylab("% of seats in the Lower Chamber") + xlab("Left-Right Dimension") +
                  labs(title = "Herrera's Pre-Electoral Coalition") +
                  mytheme2

b <- herrera_all %>%
     ggplot(aes(ideology, legseat, label = party)) +
     geom_point(aes(colour = "black", fill = "black"), size = 6, alpha = 0.85) +
     scale_color_manual(values = c("black", "black")) +
     scale_fill_manual(values = c("black", "black")) +
     geom_text(aes(x=ideology, y=legseat+0.03, label = ifelse(party == "MEP", "", party), size = 4, family = "regular")) +
     geom_text(aes(x=ideology, y=legseat-0.03, label = ifelse(party == "MEP", party, ""), size = 4, family = "regular")) +
     coord_cartesian(expand=T) +
     scale_x_continuous(limits = c(0, 11),
                        breaks = seq(1, 10, by = 1)) +
     scale_y_continuous(breaks = seq(00.10, 00.50, by = 00.10),
                        limits = c(-0.02, 00.50),
                        labels = scales::percent_format(accuracy = 1)) +
     ylab("% of seats in the Lower Chamber") + xlab("Left-Right Dimension") +
     labs(title = "Venezuelan Party System as of 1978") + 
     mytheme2
  
  grid.arrange(a, b, ncol=2)

#Balladares
#Once again, you should have the V-party dataset to reproduce the visualisation below.
vdem <- read.csv("V-Dem-CPD-Party-V2.csv", sep=",") %>%
        filter(country_name == "Panama") %>%
        filter(year == 1994) %>%
        select(v2paenname, v2pashname, v2paseatshare, v2pariglef_osp)
  
add <- tribble(
       ~v2paenname, ~v2pashname, ~v2paseatshare, ~v2pariglef_osp,
       "Partido Liberal Republicano", "LIBRE", 2.777778, 4.458,
       "Partido Laborista", "PALA", 1.388889, 4.458,
)
  
vdem1 <- vdem %>% 
         rbind(add) %>%
         mutate(ideology = 1+(19/6)*v2pariglef_osp)
  
  
balladares <- vdem1 %>%
              mutate(ideology_vd2 = ideology/2) %>%
              filter(v2pashname == "PRD" | v2pashname == "LIBRE" | v2pashname == "PALA") %>%
              mutate(v2paseatshare = v2paseatshare/100)
  
  
balladares_all  <- vdem1 %>%
                   drop_na() %>%
                   mutate(ideology_vd2 = ideology/2) %>%
                   mutate(v2paseatshare = v2paseatshare/100)

polas <- balladares_all %>%
         mutate(all = sum(v2paseatshare)) %>%
         mutate(pct = v2paseatshare/all) %>%
         mutate(mlr = pct*ideology_vd2) %>%
         mutate(general = sum(mlr)) %>%
         mutate(polar = (ideology_vd2-general)/9.5 * 10) %>%
         mutate(polar = sqrt((polar*pct)^2)) %>%
         mutate(polar = sum(polar)) %>%
         select(polar) %>% unique() %>% pull()

a <- balladares  %>% ggplot(aes(ideology_vd2, v2paseatshare, label = v2pashname)) +
     geom_point(aes(colour = "black", fill = "black"), size = 6, alpha = 0.85) +
     scale_color_manual(values = c("black", "black")) +
     scale_fill_manual(values = c("black", "black")) +
     geom_text(aes(x=ideology_vd2, y=v2paseatshare+0.03, label = ifelse(v2pashname %in% c("PALA","LIBRE"), "", v2pashname)), size = 4) +
     geom_text(aes(x=ideology_vd2+1.1, y=v2paseatshare, label = ifelse(v2pashname == "PALA", v2pashname, "")), size = 4) +
     geom_text(aes(x=ideology_vd2+1.3, y=v2paseatshare+0.01, label = ifelse(v2pashname == "LIBRE", v2pashname, "")), size = 4) +
     coord_cartesian(expand=T) +
     scale_x_continuous(limits = c(0, 11),
                        breaks = seq(1, 10, by = 1)) +
     scale_y_continuous(breaks = seq(00.10, 00.50, by = 00.10),
                        limits = c(0.00, 00.50),
                        labels = scales::percent_format(accuracy = 1)) +
     ylab("% of seats in the Only Chamber") + xlab("Left-Right Dimension") +
     labs(title = "Balladares's Pre-Electoral Coalition") +
     mytheme2
  
b <- balladares_all %>% mutate(ministerial_post = ifelse(v2pashname %in% c("PS"), 0, 1)) %>%
     ggplot(aes(ideology_vd2, v2paseatshare, label = v2pashname, colour = as.factor(ministerial_post), fill = as.factor(ministerial_post))) +
     geom_point(size = 6, alpha = 0.85) +
     scale_color_manual(values = c("#a9a9a9", "black")) +
     scale_fill_manual(values = c("#a9a9a9", "black")) +
     geom_text(aes(x=ideology_vd2, y=v2paseatshare+0.03, label = ifelse(v2pashname %in% c("PP-PDC", "LIBRE", "PALA", "MOLIRENA"), "", v2pashname)), size = 4) +
     geom_text(aes(x=ideology_vd2, y=v2paseatshare+0.02, label = ifelse(v2pashname == "PP-PDC", v2pashname, "")), size = 4) +
     geom_text(aes(x=ideology_vd2+0.5, y=v2paseatshare+0.04, label = ifelse(v2pashname == "MOLIRENA", v2pashname, "")), size = 4) +
     geom_text(aes(x=ideology_vd2+1.3, y=v2paseatshare+0.01, label = ifelse(v2pashname == "LIBRE", v2pashname, "")), size = 4) +
     geom_text(aes(x=ideology_vd2+1.1, y=v2paseatshare, label = ifelse(v2pashname == "PALA", v2pashname, "")), size = 4) +
     scale_x_continuous(limits = c(0, 11),
                        breaks = seq(1, 10, by = 1)) +
     scale_y_continuous(breaks = seq(00.10, 00.50, by = 00.10),
                        limits = c(0.00, 00.50),
                        labels = scales::percent_format(accuracy = 1)) +
     ylab("% of seats in the Only Chamber") + xlab("Left-Right Dimension") +
     labs(title = "Panamanian Party System as of 1994") + 
     mytheme2

grid.arrange(a, b, ncol=2)

#Checking on the legislative seats of excluded parties
#You should have the dataset used for Borges et al. (forthcoming) for this
nr_p <- read.csv("NewRight_partylevel_dataset_v1.2.csv", sep = ";")

nr_p1 <- nr_p %>%
         select(countrycode, partycode, legseat) %>%
         drop_na() %>%
         filter(legseat < 1.0) %>%
         group_by(countrycode) %>%
         mutate(summing = sum(legseat)) %>%
         select(countrycode, summing) %>%
         unique() %>%
         arrange(desc(summing))

bra2018 <- masters %>%
           filter(countrycode2 == "Brazil2018") %>%
           filter(composition == "PT PP PMDB PSD PR PSB REP PFL PSDB PDT PTB PCdoB NOVO PSOL PSL") %>%
           select(composition) %>%
           unique() %>%
           pull()

nr_p2 <- nr_p %>%
         filter(countrycode == "Brazil2018") %>%
         filter(party %notin% c("PT", "PP", "PMDB", "PSD", "PR", "PSB", "PRB", "PFL", "PSDB", "PDT", "PTB", "PCdoB", "NOVO", "PSOL", "PSL")) %>%
         mutate(summing = sum(legseat)) %>%
         select(summing) %>%
         unique() %>%
         pull()

#Descriptive Statistics
masters_appendix1 <- masters %>%
                     mutate(empgov_mean = mean(empgov)) %>%
                     mutate(empgov_sd = sd(empgov)) %>%
                     mutate(empgov_min = min(empgov)) %>%
                     mutate(empgov_max = max(empgov)) %>%
                     select(empgov_mean, empgov_sd, empgov_min,
                            empgov_max) %>%
                     unique()

masters_appendix2 <- masters %>%
                     mutate(minority_mean = mean(minority)) %>%
                     mutate(minority_sd = sd(minority)) %>%
                     mutate(minority_min = min(minority)) %>%
                     mutate(minority_max = max(minority)) %>%
                     select(minority_mean, minority_sd, minority_min,
                            minority_max) %>%
                     unique()

masters_appendix3 <- masters %>%
                     mutate(numpar_mean = mean(numpar)) %>%
                     mutate(numpar_sd = sd(numpar)) %>%
                     mutate(numpar_min = min(numpar)) %>%
                     mutate(numpar_max = max(numpar)) %>%
                     select(numpar_mean, numpar_sd, numpar_min,
                            numpar_max) %>%
                     unique()

masters_appendix4 <- masters %>%
                     mutate(division2_mean = mean(division2)) %>%
                     mutate(division2_sd = sd(division2)) %>%
                     mutate(division2_min = min(division2)) %>%
                     mutate(division2_max = max(division2)) %>%
                     select(division2_mean, division2_sd, division2_min,
                            division2_max) %>%
                     unique()

masters_appendix5 <- masters %>%
                     mutate(medianparty_mean = mean(medianparty)) %>%
                     mutate(medianparty_sd = sd(medianparty)) %>%
                     mutate(medianparty_min = min(medianparty)) %>%
                     mutate(medianparty_max = max(medianparty)) %>%
                     select(medianparty_mean, medianparty_sd, medianparty_min,
                            medianparty_max) %>%
                     unique()

masters_appendix6 <- masters %>%
                     mutate(extremeparties_mean = mean(extremeparties)) %>%
                     mutate(extremeparties_sd = sd(extremeparties)) %>%
                     mutate(extremeparties_min = min(extremeparties)) %>%
                     mutate(extremeparties_max = max(extremeparties)) %>%
                     select(extremeparties_mean, extremeparties_sd, extremeparties_min,
                            extremeparties_max) %>%
                     unique()

masters_appendix7 <- masters %>%
                     mutate(secfin_mean = mean(secfin)) %>%
                     mutate(secfin_sd = sd(secfin)) %>%
                     mutate(secfin_min = min(secfin)) %>%
                     mutate(secfin_max = max(secfin)) %>%
                     select(secfin_mean, secfin_sd, secfin_min,
                           secfin_max) %>%
                     unique()

masters_appendix8 <- masters %>%
                     mutate(pec_mean = mean(pec)) %>%
                     mutate(pec_sd = sd(pec)) %>%
                     mutate(pec_min = min(pec)) %>%
                     mutate(pec_max = max(pec)) %>%
                     select(pec_mean, pec_sd, pec_min,
                            pec_max) %>%
                     unique()

masters_appendix9 <- masters %>%
                     drop_na(polariz_st) %>%
                     mutate(polariz_st_mean = mean(polariz_st)) %>%
                     mutate(polariz_st_sd = sd(polariz_st)) %>%
                     mutate(polariz_st_min = min(polariz_st)) %>%
                     mutate(polariz_st_max = max(polariz_st)) %>%
                     select(polariz_st_mean, polariz_st_sd, polariz_st_min,
                            polariz_st_max) %>%
                     unique()


##Robustness Analyses----
#Congruent PECs and Enlarged PECs
masters_agg <- masters %>%
               filter(version == 2) %>%
               mutate(pecpar = case_when(countrycode2 == "Argentina2013" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%
               mutate(pecpar = case_when(countrycode2 == "Argentina2017" & composition == "UCR PRO" ~ 1,
                                        T ~ as.numeric(pecpar))) %>%
               mutate(pecpar = case_when(countrycode2 == "Argentina2019" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%     
               mutate(pecpar = case_when(countrycode2 == "El Salvador2012" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%                           
               mutate(pecpar = case_when(countrycode2 == "El Salvador2015" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%  
               mutate(pecpar = case_when(countrycode2 == "El Salvador2018" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%   
               mutate(pecpar = case_when(countrycode2 == "Honduras2013" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%     
               mutate(pecpar = case_when(countrycode2 == "Honduras2017" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%     
               mutate(pecpar = case_when(countrycode2 == "Panama2014" & composition == "PPPDC PPA" ~ 1,
                                         T ~ as.numeric(pecpar))) %>%     
               mutate(pecpar = case_when(countrycode2 == "Panama2019" & composition == "MOLIRENA PRD" ~ 1,
                                         T ~ as.numeric(pecpar))) %>%  
               mutate(pecpar = case_when(countrycode2 == "Uruguay2019" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%  
               mutate(pecpar = case_when(countrycode2 == "Argentina2015" & composition == "UCR ARI PRO" ~ 1,
                                        T ~ as.numeric(pecpar))) %>%  
               mutate(pecpar = case_when(countrycode2 == "Chile2013" & composition == "PDC PPD PCCh PRSD PSCh" ~ 1,
                                        T ~ as.numeric(pecpar))) %>%  
               mutate(pecpar = case_when(countrycode2 == "Chile2017" & composition == "UDI EP RN" ~ 1,
                                        T ~ as.numeric(pecpar))) %>%  
               mutate(pecpar = case_when(countrycode2 == "Costa Rica2014" ~ 0,
                                        T ~ as.numeric(pecpar))) %>%  
               mutate(pecpar = case_when(countrycode2 == "Costa Rica2018" ~ 0,
                                        T ~ as.numeric(pecpar))) %>%  
               mutate(pecpar = case_when(countrycode2 == "Uruguay2014" ~ 0,
                                        T ~ as.numeric(pecpar))) %>%          
               mutate(pecpar = case_when(countrycode2 == "Bolivia2014" ~ 0,
                                        T ~ as.numeric(pecpar))) %>%          
               mutate(pecpar = case_when(countrycode2 == "Colombia2018" ~ 0,
                                        T ~ as.numeric(pecpar))) %>%          
               mutate(pecpar = case_when(countrycode2 == "Colombia2014" & composition == "PLC PCR PU" ~ 1,
                                         T ~ as.numeric(pecpar))) %>%          
               mutate(pecpar = case_when(countrycode2 == "Nicaragua2011" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%          
               mutate(pecpar = case_when(countrycode2 == "Nicaragua2016" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%             
               mutate(pecpar = case_when(countrycode2 == "Brazil2014" & composition == "PMDB PSD PP PR PRB PDT PROS PCdoB PT" ~ 1,
                                         T ~ as.numeric(pecpar))) %>%             
               mutate(pecpar = case_when(countrycode2 == "Brazil2018" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%                 
               mutate(pecpar = case_when(countrycode2 == "Brazil2010" & composition == "PMDB PR PSB PDT PSC PCdoB PRB PT" ~ 1,
                                         T ~ as.numeric(pecpar))) %>%    
               mutate(pecpar = case_when(countrycode2 == "Dominican Republic2010" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%    
               mutate(pecpar = case_when(countrycode2 == "Dominican Republic2012" ~ 0,
                                         T ~ as.numeric(pecpar))) %>%    
               mutate(pecpar = case_when(countrycode2 == "Dominican Republic2016" & composition == "PRD PLD" ~ 1,
                                         T ~ as.numeric(pecpar))) %>%
               mutate(pecpar = replace_na(pecpar, 0)) %>% 
               mutate(pec = as.numeric(pec),
                      pec = ifelse(pec == 2, 1, 0), 
                      pecpar = as.numeric(pecpar)) %>%
               mutate(pecext = pec - pecpar)

master_org <- masters %>%
              filter(version == 1) %>%
              mutate(pec = as.numeric(pec),
                     pec = ifelse(pec == 2, 1, 0)) %>%
              rbind(masters_agg) %>%
              filter(majsit == 0) 

model_alt_1 <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pecpar + pecpar*polariz_st, 
  data=master_org
)


model_alt_2 <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pecext + pecext*polariz_st, 
  data=master_org
)

model_alt_3 = clogit(empgov ~ minority + numpar + division2 + medianparty +
                extremeparties + secfin + pecpar + pecext + pecpar*polariz_st + pecext*polariz_st +  
                strata(cab_id), data=master_org)

stargazer(model_alt_1, model_alt_2, model_alt_3, 
          type = "latex",
          order = c(1,2,3,4,5,6,7,9,8,11),
          keep.stat = "n", 
          covariate.labels = c("Minority", "Number of Parties", 
                               "Ideological Division", "Median Party", 
                               "Extreme Parties", "Runner-up Party",
                               "Congruent Pre-Electoral Coalition (PEC)", 
                               "Enlarged Pre-Electoral Coalition (PEC)",
                               "Congruent PEC * Legislative Polarisation",
                               "Enlarged PEC * Legislative Polarisation", "exc"),
          model.names = F,
          column.separate = c(1, 1, 1),
          dep.var.caption	= "",
          dep.var.labels.include = F,
          no.space=TRUE,
          add.lines=list(c("Number of Alternative Cabinets","147,736", "147,736", "147,736"),
                         c("Log Likelihood", "-280.698", "-303.854", "-250.693")))




#Original
first_test_official <- masters %>%
                       filter(version == 1,
                              majsit == 0)

m1 <- mclogit(
  empgov|cab_id ~ minority + uppermin + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=first_test_official
)

#Bicameral Systems Original
first_test_officialb <- masters %>%
                        filter(version == 1,
                               majsit == 0,
                               unicameral != 1)

m2 <- mclogit(
  empgov|cab_id ~ minority + uppermin + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=first_test_officialb
)

#Bicameral Systems overall
first_test_officialc <- masters %>%
  filter(majsit == 0,
         unicameral != 1)

m3 <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=first_test_officialc
)

#Baker and Greene (2011)
second_test <- masters %>%
               filter(majsit == 0)

m4 <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + pec + secfin + pec*polariz_st_BG,
  data=second_test
)

#ENPP > 2.5
third_test <- masters %>%
              filter(majsit == 0,
                     ENPP > 2.5)

m5 <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + pec + secfin + pec*polariz_st,
  data=third_test
)

#Primaries
full <- read.csv("govfor_rr.csv", sep = ",") %>%
  mutate(pec = as.factor(pec)) 

masters_wo_primary <- full %>%
                      filter(majsit == 0, 
                             primary_election == 1)

m6 = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + pec + secfin + pec*polariz_st,
  data=masters_wo_primary
)

#Proportional
masters_proportional <- masters %>%
                        filter(legislative_type == "Proportional",
                               majsit == 0)

m7 = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + pec + secfin + pec*polariz_st,
  data=masters_proportional
)

##D'Hondt
masters121 <- masters %>% select(formula) %>% unique()

masters_dhondt <- masters %>%
  filter(formula == "D<92>Hondt" | formula == "D\x92Hondt") %>%
  filter(majsit == 0)

m8 = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + pec + secfin + pec*polariz_st,
  data=masters_dhondt
)

##PSI
median(full$v2xps_party, na.rm = T)
hist(full$v2xps_party)

bas <- full %>%
       filter(v2xps_party >= 0.724) %>%
       filter(majsit == 0)


m9 <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=bas
)

##elec_volat
median(full$elec_volat, na.rm = T)
hist(full$elec_volat)

bas1 <- full %>%
       filter(elec_volat <= 2) %>%
       filter(majsit == 0)

m10 <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=bas1
)


stargazer(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10,
          type = "latex",
          order = c(1,2,3,4,5,6,7,9,10,8),
          covariate.labels = c("Lower Chamber Minority", "Upper Chamber Minority",
                               "Number of Parties", "Ideological Division", "Median Party",
                               "Extreme Parties", "Runner-up Party", "Pre-Electoral Coalition (PEC)",
                               "PEC * Legislative Polarisation", "PEC * Legislative Polarisation (BG)"),
          dep.var.labels.include = F)

#Clientelism
#The dataset is available at https://sites.duke.edu/democracylinkage/data/.
master_countries <- masters %>% 
                    select(country_name_long) %>% 
                    mutate(country_name_long = ifelse(country_name_long == "Dominican Republic", "Dom. Rep.", country_name_long)) %>%
                    unique() %>% 
                    pull()

clientelism <- read_dta("countrylevel_20130907.dta")

client_avg <- clientelism %>% select(b6) %>% mutate(median = median(b6)) %>% select(median) %>% unique() %>% pull()

client1 <- clientelism %>%
           filter(country %in% master_countries) %>%
           select(country_name_long = country, b6)

masters_client <- masters %>%
                  mutate(country_name_long = ifelse(country_name_long == "Dominican Republic", "Dom. Rep.", country_name_long)) %>%
                  left_join(client1, by = "country_name_long") %>%
                  filter(b6 >= client_avg) %>%
                  filter(majsit == 0) 


m11 <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_client
)

#programmaticness
#Available at https://v-dem.net/data/dataset-archive/. Version 13 used below.
vdem <- read.csv("V-Dem-CY-Full+Others-v13.csv", sep = ",") 

vdem_p <- vdem %>%
          select(country_name, year, v2psprlnks_osp) %>%
          mutate(countrycode2 = paste0(country_name, year)) %>%
          select(countrycode2, v2psprlnks_osp)

masters_program <- masters %>%
                   left_join(vdem_p, by = "countrycode2") 

median(masters_program$v2psprlnks_osp, na.rm = T)
hist(masters_program$v2psprlnks_osp)


masters_program_program <- masters_program %>%
                           filter(v2psprlnks_osp >= 2.967) %>%
                           filter(majsit == 0) 

m12 <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_program_program
)

masters_program_cleint <- masters_program %>%
                           filter(v2psprlnks_osp < 2.967) %>%
                           filter(majsit == 0) 

m13 <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_program_cleint
)


stargazer(m11, m12, m13, 
          type = "latex",
          order = c(1,2,3,4,5,6,7,8), 
          covariate.labels = c("Minority", "Number of Parties", "Ideological Division", "Median Party",
                               "Extreme Parties", "Runner-up Party", "Pre-Electoral Coalition (PEC)",
                               "PEC * Legislative Polarisation"),
          column.labels   = c("DALP", "Programmatic - V-Dem" , "Clientelistic - V-Dem"),
          #model.names = F,
          dep.var.labels.include = F,
          dep.var.caption	= "",
          no.space=TRUE,
          add.lines=list(c("Number of Alternative Cabinets","145,672", "106,320", "41,416"),
                         c("Log Likelihood", "-205.265", "-72.864", "-154.700")))


#liberal-conservative dimension
#Dataset for Borges et al. (forthcoming)
nr <- read.csv("NewRight_systemlevel_dataset_v1.21.csv", sep=";") 

nr1 <- nr %>%
  select(countrycode2, polz_libcons_st)

m1 <- masters %>%
  left_join(nr1, by = "countrycode2") %>%
  filter(majsit == 0)

model_lib_con <- mclogit(empgov|cab_id ~ minority + numpar + division2 + medianparty +
                           extremeparties + secfin + pec + pec*polz_libcons_st,
                         data=m1)

stargazer(model_lib_con,
          type = "latex",
          order = c(1,2,3,4,5,6,7,8),
          covariate.labels = c("Lower Chamber Minority", 
                               "Number of Parties", "Ideological Division", "Median Party",
                               "Extreme Parties", "Runner-up Party", "Pre-Electoral Coalition (PEC)",
                               "PEC * Legislative Polarisation (Liberal-Conservative Dimension)"),
          model.names = F,
          dep.var.caption	= "",
          #dep.var.labels.include = F,
          no.space=TRUE,
          add.lines=list(c("Number of Alternative Cabinets","143,164"),
                         c("Log Likelihood", "-216.647")))

#ideological distance between presidential and median parties
pp_mp <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st + pec*polariz_st*diff_pp_mp,
  data=freude
)

stargazer(pp_mp,
          type = "latex",
          order = c(1,2,3,4,5,6,7,8),
          covariate.labels = c("Lower Chamber Minority", 
                               "Number of Parties", "Ideological Division", "Median Party",
                               "Extreme Parties", "Runner-up Party", "Pre-Electoral Coalition (PEC)",
                               "PEC * Legislative Polarisation", "PEC * Diff. PP - MP", "PEC * Leg. Polar. * Diff. PP - MP"),
          model.names = F,
          column.labels   = c("Model 1"),
          dep.var.labels.include = F,
          no.space=TRUE,
          add.lines=list(c("Number of Alternative Cabinets","139,032"),
                         c("Log Likelihood", "-234.455")))

#presidential powers
prespow <- mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st + pec*polariz_st*prespow1,
  data=freude)

stargazer(prespow,
          type = "latex",
          order = c(1,2,3,4,5,6,7,8),
          covariate.labels = c("Minority", 
                               "Number of Parties", "Ideological Division", "Median Party",
                               "Extreme Parties", "Runner-up Party", "Pre-Electoral Coalition (PEC)",
                               "PEC * Legislative Polarisation", "PEC * Presidential Power", "PEC * Leg. Polar. * Pres. Power"),
          model.names = F,
          column.labels   = c("Model 1"),
          dep.var.labels.include = F,
          no.space=TRUE,
          add.lines=list(c("Number of Alternative Cabinets","147,736"),
                         c("Log Likelihood", "-243.804")))

ent <- plot_cme(prespow, variables = "pec", condition = c("polariz_st", "prespow1"), type = "link", draw = F) 

plot_cme(prespow, variables = "pec", condition = c("polariz_st", "prespow1"), type = "link") +
  ggplot2::geom_line(data = ent, ggplot2::aes(x = polariz_st, 
                                              y = estimate,
                                              color = factor(prespow1)), size = 0.8) +
  ggplot2::geom_ribbon(data = ent %>% filter(prespow1 == 0.257), alpha = 0.5, ggplot2::aes(x = polariz_st, 
                                                                                           y = estimate, 
                                                                                           ymin = conf.low, 
                                                                                           ymax = conf.high,
                                                                                           fill = prespow1)) +
  ggplot2::geom_ribbon(data = ent %>% filter(prespow1 == 0.452), alpha = 0.5,  ggplot2::aes(x = polariz_st, 
                                                                                            y = estimate, 
                                                                                            ymin = conf.low, 
                                                                                            ymax = conf.high,
                                                                                            fill = prespow1)) +
  ggplot2::geom_ribbon(data = ent %>% filter(prespow1 == 0.486), alpha = 0.2, ggplot2::aes(x = polariz_st, 
                                                                                           y = estimate, 
                                                                                           ymin = conf.low, 
                                                                                           ymax = conf.high,
                                                                                           fill = prespow1)) +
  ggplot2::geom_ribbon(data = ent %>% filter(prespow1 == 0.57), alpha = 0.2, ggplot2::aes(x = polariz_st, 
                                                                                          y = estimate, 
                                                                                          ymin = conf.low, 
                                                                                          ymax = conf.high,
                                                                                          fill = prespow1)) +
  geom_hline(yintercept=0, linetype="dashed", color = "red", size=2) +
  xlab("Legislative Polarisation") + ylab("Marginal Effect of Pre-Electoral Coalitions on Cabinet Formation") + guides(fill = guide_legend(title="Presidential \nPower"),
                                                                                                                       colour = guide_legend(title="Presidential \nPower")) +
  
  theme(
    ##Background
    panel.border = element_rect(color = "black", fill = NA),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "#ffffff", linetype = 'blank'),
    panel.background = element_rect(fill = "#ffffff"),
    ##Axes
    axis.title = element_text(size=15, family="regular", color="black"),
    axis.text.y = element_text(size=10, family="regular", color="black"),
    axis.text.x = element_text(size=10, family="regular", color="black"))


#Iterative Exclusions
masters <- masters %>% filter(majsit == 0)

masters_wo_arg <- masters %>% filter(country_name_long != "Argentina")  

arg = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_wo_arg
)

masters_wo_bra <- masters %>% filter(country_name_long != "Brazil")  

bra = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + pec + secfin + pec*polariz_st,
  data=masters_wo_bra
)

masters_wo_chi <- masters %>% filter(country_name_long != "Chile")  

chi = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + pec + secfin + pec*polariz_st,
  data=masters_wo_chi
)

masters_wo_cos <- masters %>% filter(country_name_long != "Costa Rica")  

cos = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + pec + secfin + pec*polariz_st,
  data=masters_wo_cos
)

masters_wo_hnd <- masters %>% filter(country_name_long != "Honduras")  

hnd = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_wo_hnd
)

masters_wo_slv <- masters %>% filter(country_name_long != "El Salvador")  

slv = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_wo_slv
)

masters_wo_ury <- masters %>% filter(country_name_long != "Uruguay")  

ury = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_wo_ury
)

masters_wo_vnz <- masters %>% filter(country_name_long != "Venezuela")  

vnz = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_wo_vnz
)

masters_wo_pnm <- masters %>% filter(country_name_long != "Panama")  

pnm = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_wo_pnm
)

masters_wo_bol <- masters %>% filter(country_name_long != "Bolivia")  

bol = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + pec + secfin + pec*polariz_st,
  data=masters_wo_bol
)

masters_wo_col <- masters %>% filter(country_name_long != "Colombia")  

col = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + pec + secfin + pec*polariz_st,
  data=masters_wo_col
)

masters_wo_nic <- masters %>% filter(country_name_long != "Nicaragua")  

nic = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_wo_nic
)

masters_wo_dom <- masters %>% filter(country_name_long != "Dominican Republic")  

dom = mclogit(
  empgov|cab_id ~ minority + numpar + division2 + medianparty +
    extremeparties + secfin + pec + pec*polariz_st,
  data=masters_wo_dom
)

stargazer(arg, bol, bra, chi, col, cos, dom, 
          type = "latex",
          covariate.labels = c("Minority", "Number of Parties", 
                               "Ideological Division", "Median Party", 
                               "Extreme Parties", "Runner-up Party",
                               "Pre-Electoral Coalition (PEC)", " PEC * Legislative Polarisation"),
          dep.var.labels.include = F, 
          dep.var.caption = c("Minority Presidents"))

model_3= clogit(empgov ~ minority + numpar + division2 + medianparty +
                  extremeparties + secfin + pec + pec*polariz_st +
                  strata(cab_id), data=masters_wo_dom)


stargazer(slv, hnd, nic, pnm, ury, vnz,  
          type = "latex",
          covariate.labels = c("Minority", "Number of Parties", 
                               "Ideological Division", "Median Party", 
                               "Extreme Parties", "Runner-up Party",
                               "Pre-Electoral Coalition (PEC)", " PEC * Legislative Polarisation"),
          dep.var.labels.include = F, 
          dep.var.caption = c("Minority Presidents"))

